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We examine the influence of the main approximations employed in density-functional theory de- 
scriptions of the solid phase of molecular hydrogen near dissociation. We consider the importance of 
nuclear quantum effects on equilibrium properties and find that they strongly influence intramolec- 
ular properties, such as bond fluctuations and stability. We demonstrate that the combination of 
both thermal and quantum effects make a drastic change to the predicted optical properties of 
the molecular solid, suggesting a limited value to dynamical, e.g., finite-temperature, predictions 
based on classical ions and static crystals. We also consider the influence of the chosen exchange- 
correlation density functional on the predicted properties of hydrogen, in particular, the pressure 
dependence of the band gap and the zero-point energy. Finally, we use our simulations to make an 
assessment of the accuracy of typically employed approximations to the calculation of the Gibbs 
free energy of the solid, namely the quasi-harmonic approximation for solids. We find that, while 
the approximation is capable of producing free energies with an accuracy of ~ 10 meV, this is not 
enough to make reliable predictions of the phase diagram of hydrogen from first-principles due to 
the small free energy differences seen between several potential structures for the solid; direct free 
energy calculations for quantum protons are required in order to make definite predictions. 

PACS numbers: 67.80.ff,63.20.dk,62.50.-p,64.70.kt 



I. INTRODUCTION 

Hydrogen at high pressure, being the most abundant 
element in the Universe, plays a prominent role in many 
scientific fields, notably in modeling of the giant planets. 
Its simplicity, i.e., being comprised of a single proton 
and electron per atom, makes it a unique system from 
a theoretical and computational standpoint. Thus, con- 
siderable attention has been devoted to understanding 
hydrogen both experimentally and theoretically, as has 
been recently reviewed^. 

Numerous exciting experimental breakthroughs have 
recently been made in the low-temperature solid phase 
of molecular hydrogen, including a controversial obser- 
vation of metallizatiorP as well as the discovery of a 
new phase, Phase rV™. Despite these advances, many 
open questions remain, such as whether metallization has 
actually been achievecP and the structure of this new 
phased Co mp utational predictions have accompanied 
these results^. However, the agreement of these pre- 
dictions with experiment is not perfect, and have not 
resolved these questions. 

To help resolve these discrepancies, a careful ex- 
amine of the approximations made in such simula- 
tions is neededd. For example, in almost all studies 
of hydrogen reported to date, especially in the solid 
phase, many properties have been calculated for perfect 
latticed^', with a f inite- temperature description based 
on classic al proto ns^^, or within a quasi-harmonic 
calculatio n 1 5 1 6 1 12 1 13 - 1 . While such procedures would be jus- 
tified if both nuclear quantum effects (NQEs) and ther- 



mal fluctuations were small, such is not the case for hy- 
drogen at the high pressures under consideration. Notice 
that even at temperatures below T — 100 K, the kinetic 
energy of the protons in the crystal is on the order ^1000 
K, because of proton zero point energy (ZPE). While 
some progress has been made in the direct treatment 
of NQEs in hydrogen from first principles recentljEHUD 
there is still considerable work to be done. 

Another major approximation arises in density- 
functional theory (DFT) studies of hydrogen, namely the 
choice of the approximate exchange-correlation density 
functional (DF). Local or semi- local DFs, such as that of 
Perdew-Burke-Ernzerhof (PBE^l have been the stan- 
dard choice in DFT simulations over the last decade. 
However, it has recently been demonstrated, at least at 
higher temperature in the liquid phase^, that the use 
of nonlocal functionals, such as that by Heyd-Scuseria- 
Ernzerhof (HSE^l which contains a fraction of exact ex- 
change or one with an improved description of dispersion 
interactions, such as vdW-DF2, significantly improve the 
description of the dissociation process. Since the neglect 
of NQEs and the deficiencies of the PBE DF partially 
compensate each other in many situations^, most cal- 
culations to date have made the assumption that this 
cancellation is accurate. As we show below, accurate 
and predictive simulations of hydrogen, in particular in 
the region where molecular dissociation and metallization 
occur, require a more rigorous treatment of electronic 
exchange and correlation effects (i.e., beyond those pro- 
vided by typical semi-local functionals like PBE). 

The main purpose of this article is to examine the effect 
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of these approximations in first-principles simulations of 
crystalline molecular hydrogen, in order to improve the 
accuracy of the method. We focus our study on the in- 
fluence of NQEs and the choice of DF on the orienta- 
tional order in the crystal and on bandgaps. We also 
make a critical assessment of some common approxima- 
tions, such as the quasi-harmonic approximation applied 
to the prediction of the relative energies of competing 
structures. We begin in Section [TT] by providing the com- 
putational details for the results that follow in Section 
TIT| Section [W] provides a brief discussion of the results, 
and ScctionTvlconcludes. 



II. COMPUTATIONAL DETAILS 

We performed first-principles simulations of hydrogen 
using DFT. Three DFs were considered: the semi-local 
PBE DF23, the HSE DF^, containing^ fraction of exact- 
exchange, and the vdW-DF2 DF22H22 capable of treat- 
ing van der Waals (vdW) interactions. Simulations using 
HSE were performed with a modified version of VASF 23 , 
while those with PBE and vdW-DF2 were performed 
with a modified version of Quantum ESPRESSO (QE^P. 
A Troullier-Martins norm-conserving pseudopotential 25 
with a core radius of r c — 0.5 a. u. was used to replace the 
bare Coulomb-potential of hydrogen in the QE simula- 
tions and a PAW^H was used in VASP. Planewave cutoffs 
of 1224 eV and 350 eV were used in these simulations, 
respectively. 

Path-integral molecular dynamics (P1MD) simulations 
were employed via the accelerated method of Ceriotti, et 
a?P3, based on a generalized Langevin dynamics (GLE) 
and the Born-Oppenheimer (BO) approximation, which 
we indicate as PI+GLE. The use of the PI+GLE method 
was carefully tested under the relevant pressure and 
temperature conditions, in order to guarantee proper 
convergence^. A time-step of 8 (a. u.) _1 was used in 
all simulations, and the Pis were discretized with a Trot- 
ter time-step no larger than 0.0003125 K _1 . After an 
equilibration period of ~ 0.25 ps, statistics were gath- 
ered for simulation times of ~1. 25-1. 75 ps, correspond- 
ing to ~ 6500-9000 time steps. A 2 3 k-point grid was 
used in simulations with the PBE and vdW-DF2, while 
a 1 x 2 x 2 grid was used in simulations with HSEP^. 
Finite-temperature effects on the electrons were taken 
into account using Fermi-Dirac smearing^. 

We performed both classical nuclei and quantum sim- 
ulations for all DFs of three of the primary candidate 
structures for the high-pressure molecular phases (in the 
region of ~300 GPap? (72c, Cmca-12, and Pbcn. All 
calculations were performed in the NVT ensemble, where 
N is the number of particles and V is the volume. Sim- 
ulation cells contained 144 atoms at temperatures be- 
tween 200 K and 500 K with pressures from 200-550 
GPa. Note that below, we use the term BOMD to refer 
to simulations that treat the protons as classical parti- 
cles, while PIMD refers to those including a full path 




FIG. 1: (Color online) The proton-proton PCFs for the C2c 
(left), Cmca-12 (center), and Pbcn (right) structures of hy- 
drogen for BOMD (dashed blue) and PIMD (solid red) sim- 
ulations, using the PBE DF. From top to bottom the PCFs 
correspond to pressures of p ~350, 300, 250, and 200 GPa. 
All simulations were performed at T = 200 K. 



integral treatment of the protons (except for quantum 
statistics of the protons). Bandgaps, at finite tempera- 
ture, were calculated by performing excited-state calcu- 
lations on 15 statistically-independent proton configura- 
tions sampled from the trajectories. Since semi-local DFs 
are well known to underestimate the bandgarj^, unless 
otherwise stated, the HSE functional was used to cal- 
culate bandgaps. In other words, the trajectories and 
optical properties were not necessarily calculated using 
the same DF, and for simplicity, in our descriptions, the 
label refers to the DF used to generate the trajectories. 



III. RESULTS 
A. Orientational-order 

We begin by considering perhaps the property of prime 
importance, the crystal structure of the solid at finite 
temperature. In particular we look at its orientational 
order. Since most studies to date have employed PBE 
in DFT studies^, we use this choice initially to examine 
NQEs, but later we will assess the influence of other DFs. 

Figure [l] shows comparisons of the pair correlation 
functions (PCFs) of hydrogen computed with BOMD 
and PIMD in the C2c, Cmca-12, and Pbcn phases. A 
marked disagreement between classical and quantum re- 
sults is clearly seen; while classical simulations produce 
PCFs with considerable structure, including the exis- 
tence of molecules with different bond lengths (as pre- 
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FIG. 2: (Color online) Lindemann ratio (top) and orienta- 
tional order parameter (bottom) of various structures of solid 
molecular hydrogen from BOMD (solid symbols) and PIMD 
(open symbols) simulations, using the PBE DF. All the sim- 
ulations were performed at a temperature of T = 200 K. 
Squares correspond to C2c, circles to Cmca-12, and triangles 
to Pbcn. 



viously reportec ! 5 * 6 * 30 * ^^ ; the quantum simulations are 
seen to have much less structure. 

Another interesting feature of Fig. [T] is the overlap be- 
tween neighboring molecules when NQEs are included. 
In the case of classical simulations, there is a clear separa- 
tion between close molecules with almost no overlap with 
the closest shell. With the inclusion of NQEs though, 
not only is the height of the molecular peak dramati- 
cally reduced, but the first minimum disappears above 
pressures of ^300 GPa. While this might suggest the 
dissociation of molecules (i.e., since there is no clear 
separation between the intramolecular distance and the 
nearest-neighbor separation), a further examination of 
the trajectories suggests stable molecules at all pressures 
considered. 

Figure [2] shows the Lindemann ratio of the molecular 
center of mass and the orientational order parameter for 
both BOMD and PIMD simulations. There are several 
ways to measure orientational order^ we measure the 
deviation of a molecule Oj during the simulation with 
respect to the perfect static lattice orientation e: 



O = 



1 N 



(1) 



(O) = 1 in the static lattice, while a solid devoid of orien- 
tational order (e. g. Phase-I of hydrogen^) gives (O) = 0. 
While the Lindemann ratio depends on structure, as ex- 
pected, it is fairly insensitive to pressure in the studied 
range: the molecules are stable and the crystals do not 
melt. Note that (O) is also insensitive to pressure in this 
range. NQEs are seen to have a moderate effect in the 
Lindemann ratio of the molecular centers, which suggests 
that the marked differences observed in the PCFs come 
largely from strong NQEs on bond fluctuations. 

Notice that there is a significant difference in the 
amount of orientational order between the Cmca-12 
structure and the other two considered when NQEs are 
included. Molecules in the Cmca-12 structure largely re- 
tain the crystalline orientation even at the temperature 
investigated. In this case, NQEs are quite consistent, 
reducing the order parameter by roughly 20%. In the 
other two structures though, the amount of order of the 
perfect crystal is much more reduced by temperature al- 
ready, while NQEs play less of a role (a further reduc- 
tion of only ~5-10%). This result suggests that NQEs 
influence structures in different ways and with different 
magnitudes. 

Having established the importance of NQEs, it is also 
important to consider the influence of DF. Figure [3] 
shows a comparison of the PCFs between PBE, HSE, and 
vdW-DF2, all obtained from PIMD simulations. Several 




where Pi is the Legendre polynomial. One obtains 



FIG. 3: (Color online) PCFs for the C2c (left), Cmca-12 
(center), and Pbcn (right) structures of hydrogen from PIMD 
simulations using PBE (solid red), HSE (dotted blue), and 
vdW-DF2 (dashed-dotted black). From top to bottom, the 
PCFs correspond to pressures of approximately p ~350, 300, 
250, and 200 GPa. All simulations were performed at T — 200 
K. 

prominent features are clear. On the one hand, vdW-DF2 
simulations result in structures with more pronounced 
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FIG. 4: (Color online) Lindemann ratio (top) and orienta- 
tional order parameter (bottom) for solid molecular hydrogen 
from PIMD simulations using vdW-DF2. Note that all the 
simulations were performed at T = 200 K. Squares corre- 
spond to C2c, circles to Cmca-12, and triangles to Pbcn. 
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FIG. 5: (Color online) Bandgap of C2c using PBE. Red 
squares are the bandgaps for the static crystal, green circles 
are band gaps for classical nuclei at 200 K, and blue triangles 
are band gaps for quantum nuclei at 200 K. Cyan crosses are 
recent experimental measurements from Goncharov, et aZ.p^. 
Only one set of experimental results is shown, in order to 
provide a scale that allows a clear comparison between the 
BOMD and PIMD simulations. See Figure[7|for a more com- 
plete set of experimental results. 



molecular features, while those with PBE lead to the 
structures with less molecular features. Although, at 
high pressures, the molecular peak becomes very weak, 
molecules do not show any tendency to dissociate. Simi- 
lar to the results shown in Fig. [2] this can be seen in the 
Lindemann ratio of the molecules and the orientational 
order parameter, shown in Fig.[4]for the simulations with 
vdW-DF2. A somewhat large, but stable, Lindemann ra- 
tio can be seen for Pbcn at the two lowest pressures, with 
similar behavior for C2c at the higher pressures. This 
suggests a distortion in these structures at finite temper- 
ature, however, the molecules stay intact. 

Note also that the orientational order parameter goes 
to zero at pressures below ^200 GPa in both the Pbcn 
and C2c structures, suggesting a possible transition to 
Phase-I somewhere between 200 and 250 GPa. A precise 
prediction of the location of the I-III phase boundary 
requires the treatment of bosonic exchange, which we are 
not considering here. 



B. Bandgaps 

Having established the influence of NQEs and DFs on 
the finite-temperature structure, we now turn to their 
influence on the bandgap. Figure [5] shows a comparison 
of the pressure dependence of the electronic bandgap for 



seen to be dramatically affected by NQEs, not only in its 
magnitude, but also in its pressure dependence. While 
calculations on static crystals result in bandgaps in rea- 
sonable agreement with experiment, the proper inclusion 
of NQEs leads to conditions where the bandgap closes at 
pressures as low as p ~200 GPa. Furthermore, for PIMD 
simulations using PBE, the bandgaps of all three struc- 
tures close below 250 GPa. Note that these res ults are 
in disagreement with experimental measurements^^lEZ] 

While we have already established a strong influence 
of the DF, the same behavior is observed in all structures 
considered and regardless of the DF. Figs. [6] and [7] show 
a comparison of the pressure dependence of the bandgap 
for PIMD simulations performed with HSE and vdW- 
DF2; note that results using PBE are not shown, since, 
as just discussed, the bandgap goes to zero between 200 
- 250 GPa in all structures. Some experimental results 
are also shown in these figures. Several points are worth 
emphasizing. On the one hand, the use of HSE does 
not lead to a significant improvement compared to PBE, 
since, while the bandgaps are finite, they are still consid- 
erably below the experimental ones. On the other hand, 
results obtained with vdW-DF2 agree very well, both the 
magnitude of the bandgap and its pressure dependence. 
These results suggest that C2c or Pbcn are more suit- 
able structures in the investigated pressure range and at 
T=200K. 
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FIG. 6: (Color online) Bandgap as a function of pressure. 
Filled symbols correspond to PIMD simulations with HSE, 
while empty symbols correspond to experimental results; lines 
represent guides to the eye. Red squares, green circles, and 
blue triangles are theoretical results for the C2c, Cmca-12 
and Pbcn structures, respectively. The experimental results 
correspond to: orange pentagrams, Loubeyre, et alW\ brown 
diamonds, Zha et a/.,, cyan asterisks, Goncharov, et aZ.p^, 
and downward black triangles, Howie, et aZ.pl 



FIG. 8: (Color online) Internal energy of Cmca-12, as a func- 
tion of temperature, at a density of r a w 1.38. Red squares 
correspond to PIMD simulations with vdW-DF2, the solid 
blue-curve corresponds to QHA results (also using vdW-DF2) 
and the dashed black-curve corresponds to the energy of an 
isolated quantum rotor, displaced to match the PIMD energy 
at 200 K and with a rotational constant of an isolated hydro- 
gen molecule of Quot ~ 87 K. 
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FIG. 7: (Color online) Bandgap as a function of pressure. 
Filled symbols correspond to PIMD simulations with vdW- 
DF2. See the caption of Figure [6] for additional details. 

C. Energetics 

While the previous discussion shows the dramatic in- 
fluence of NQEs on the structural and optical properties 
of the solid, it is important to emphasize that no attempt 
was made to correct our BOMD data for quantum effects, 
as is typically done for other properties such as pressures 
and free energies. The common approach, especially in 
the study of high-pressure hydrogen, is to use the quasi- 
harmonic approximation (QHA) for solids 38 to incorpo- 
rate ZPE corrections to an otherwise classical treatment 
of the nuclei. While this method is fairly accurate in 
many materials, leading to reasonable predictions when 
compared to experiment!^, its applicability to light el- 



ements at high pressure is problematical. In the case 
of solid molecular hydrogen, the combination of orienta- 
tional order, large amplitude fluctuations, and large an- 
harmonic effects (both classical and quantum) make the 
application of this approach questionable, when taking 
into account the accuracy required for the correct deter- 
mination of crystal structures^. 

Figure [H] shows a comparison of the temperature de- 
pendence of the total energy of Cmca-12 at ^265 GPa 
between PIMD results using vdW-DF2 and the QHA. A 
clear discrepancy in the overall temperature dependence 
is apparent. While the PIMD results show a clear linear 
dependence in the regime studied, the QHA results show 
a very small effect up to ~200 K. The linear temperature- 
dependence of the energy in this regime can be under- 
stood by considering the fact that a hydrogen molecule 
has a rotational constant of approximately Ofl t ~ 87 
K. This implies that the rotational component to the 
heat capacity acquires its classical value of Cy ot = ks 
above ~100K in the case of distinguishable nuclei con- 
sidered here^. Figure [8] also shows the energy of an 
isolated quantum rotor (again assuming distinguishable 
particles), displaced to match the PIMD energy at 200 
K. Notice the excellent agreement in the resulting tem- 
perature dependence, supporting the fact that molecular 
rotations behave classically in this regime. On the other 
hand, both intra-molecular and center-of-mass vibrations 
show strong quantum behavior in this regime. The re- 
sults obtained with the QHA, on the other hand, do not 
capture this behavior, resulting in an incorrect temper- 
ature dependence of the energy. These results put into 
question the recent predictions of the transition lines be- 
tween Phases III and IV in molecular hydrogen based on 
QHA calculations^. 
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FIG. 9: (Color online) Comparison between PIMD and the 
QHA of the thermal and quantum contribution to the internal 
energy of hydrogen in C2c and Cmca-12 at K (top) and 200 
K (bottom). Symbols represent results with PIMD and lines 
the QHA. Note that in both cases, PBE was used. 
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FIG. 10: (Color online) Comparison between PIMD and the 
the QHA of the thermal and quantum contribution to the 
pressure component of the enthalpy in C2c and Cmca-12. 
Symbols represent results with PIMD and lines with the QHA. 
Note that in both cases, PBE was used. 



Figures [9] and [TU] show a comparison between PIMD 
and the QHA of the combined thermal and quantum con- 
tributions to the internal energy and to the pressure com- 
ponent of the enthalpy (PV) using PBE. Only results 
for C2c and Cmca-12 are considered, since the ground 
state structure of Pbcn displays imaginary phonons at 



T = K. For any thermodynamic quantity A, we de- 
fine its thermal and quantum contribution as A Ax (T) — 
A X (T) - A crystal , X represents either PIMD or QHA, 
and A crys tai is the property calculated at K. Note also 
that the PIMD results at K are an estimation of the 
ZPE using the quantum rotor model, as shown in Fig. [8] 
While, as it has been demonstrated above, PBE is not 
expected to provide a good description of solid molecular 
hydrogen close to the metallization, these results provide 
a useful benchmark to measure the expected accuracy 
of the QHA. Unfortunately though, there is a consider- 
able discrepancy between PIMD and the QHA, with the 
latter producing a large overestimate of the zero-point 
contribution to both energies and enthalpies. 



Figure 1 1 shows a comparison of the thermal and quan- 
tum contribution to the internal energy of the structures 
from PIMD simulations, as a function of volume at 200 
K. Note that, in addition, the error of the results if the 
QHA were to be used is presented in the inset. Notice 
that the magnitude of the energy is dependent on DF, 
with PBE producing the smallest contribution and HSE 
the largest. Also notice that the QHA error is dependent 
on both structure and functional, which eliminates the 
possibility of combining results (e.g., ground-state ener- 
gies and ZPE estimates) with different DFs to reduce 
the expense of the computations. In fact, the difference 
in energies between DFs is considerably larger than that 
between structures, for any given DF. 

While free-energy calculations have not been at- 
tempted in this work, the results presented nonetheless 
allow for a comparison of the enthalpy of solid at finite 
temperature, and also to assess the relative accuracy of 



the QHA with vdW-DF2. Figure 12 shows the enthalpy 
of all three structures calculated using PIMD with vdW- 
DF2, along with results using the QHA. Enthalpies are 
plotted relative to a polynomial fit to that of C2c. From 
these calculations, one cannot distinguish between the 
C2c and the Pbcn structure, since their enthalpies are 
within error bars of each other; also, one needs the rela- 
tive entropies to determine the stable structure. In addi- 
tion, Cmca-12 has a higher enthalpy over the considered 
temperature and density range. The relative accuracy 
of the QHA on the enthalpy is only around 5-10 meV. 
While this is sufficient to establish a small list of can- 
didate structures, it is not enough to make a definite 
prediction of the phase diagram. 

According to figure|4j the Cmca-12 structure has lower 
Lindemann ratio and higher orientational order than the 
other two structures, a fact which signals a lower en- 
tropy for Cmca-12. Although a precise statement about 
the relative stability of the considered structures can only 
be based on a free-energy calculation, a lower entropy to- 
gether with a higher enthalpy (see figure 12 ) suggests the 
Cmca-12 structure to be less favorable than the others 
over the entire pressure range considered. 
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FIG. 11: (Color online) Thermal and quantum contribution to 
the energy of hydrogen from PIMD simulations. From top to 
bottom the results correspond to: HSE, vdW-DF2, and PBE. 
Red squares, black circles, and blue triangles correspond to 
results for C2c, Cmca-12, and Pbcn, respectively. Insets show 
the error of the QHA, if it were to be employed. 



IV. DISCUSSION 

The results above clearly show that in order to produce 
an accurate first-principles description of solid molecu- 
lar hydrogen, especially close to metallization, not only 
does one need to take into account the quantum nature 
of the protons, but one also needs to go beyond stan- 
dard semi-local descriptions of the electronic structure in 
DFT. The large ZPE of the protons increases the mag- 
nitude of molecular vibrations to the point that many 
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FIG. 12: (Color online) Enthalpy of several structures of 
hydrogen at 200 K, relative to a polynomial fit to that of 
C2c. Symbols represent PIMD simulations with vdW-DF2 
and lines results from the QHA. 



structural features that are prominent in classical treat- 
ments end up being "washed out" . This can be under- 
stood by considering that the proper inclusion of NQEs 
leads to nuclear kinetic energies of ^1000 K, even at tem- 
peratures of only 200 K or lower, which are not taken 
into account in classical simulations. In other words, the 
classical picture of a molecular bond at low temperatures 
is very different from the quantum description; this can 
lead to artificially low displacements in such simulations. 
While the magnitude of the ZPE can be estimated from 
the QHA, its effects on the dynamical properties of the 
solid and its stability are not easy to estimate accurately 
with perturbative methods. 

The results also show the failure of PBE in the cor- 
rect description of molecular hydrogen close to disso- 
ciation. While it has been known for some time that 
PBE underestimates the bandgap in excited state cal- 
culations, there remains a wide-spread expectation that 
the nuclear distribution produced by this DF should be 
nonetheless accurate. It is also believed that calculations 
with PBE should produce an accurate description of the 
optical properties of hydrogen, as long as a more accu- 
rate method, like hybrid functionals or GW, is used to 
calculate the gap. This has, in fact, been u sed to predict 
excellent agreement in calculated bandgap d 37 * 43 l This is 
shown in figure 5 where a good agreement for the band 
gap of a perfect crystal predicted by PBE-DF AIRSS is 
obtained. However, this agreement is accidental as we 
have shown; even the use of HSE to sample the ionic dis- 
tribution leads to poor results. Only by incorporating 
non-local correlation capable of describing dispersion in- 
teractions provides a reasonable description of the optical 
properties of the solid obtained. Note that this is consis- 
tent with recent calculations performed in the liquid^. 
Unfortunately, the bandgaps of both C2c and Pbcn were 
very similar, which does not allow us to use these results 
to suggest any insight into the actual structure of Phase 
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III or IV. This is a topic of future work. 

Finally, the need to go beyond the QHA to estimate 
the impact of NQEs was demonstrated. While the QHA 
significantly reduces the magnitude of the errors com- 
pared to a purely classical description of the protons, 
its accuracy is probably not enough to provide predic- 
tive results in some cases, e.g., the phase diagram. The 
relative energy differences between competing structures 
using the QHA is very small, with up to half a dozen 
structures separated by energies below 10 meV 12 . While 
this approximation is certainly useful in identifying po- 
tential candidate structures, our results show that it is 
incapable of producing results with an accuracy of ~10 
meV, in some cases the error is larger. 

Looking forward, since most of the work done to date 
in the field of high-pressure hydrogen has employed PBE 
and the results above suggest that vdW-DF2 produces 
the best description of solid molecular hydrogen among 
the DFs considered in this work, it is important not only 
to revisit the problem of structure prediction using vdW- 
DF2, but it is also important to rigorously assess its accu- 
racy by comparison to more accu rate m any-body meth- 
ods, such as quantum Monte CarlcW^HUl j n addition, a 
careful study that includes accurate free-energy methods 
for quantum nuclei is needed in order to make accurate 
predictions of the correct structure ordering at finite tem- 
perature. 



V. SUMMARY AND CONCLUSIONS 

In summary, the results shown in this article clearly il- 
lustrate two important limitations of current practices in 



the theoretical study of hydrogen: on the one hand, stan- 
dard semi-local DFs like PBE lead to a poor description 
of high-pressure hydrogen^ and, on the other hand, the 
neglect of NQEs leads to predictions with very limited 
validity. While it has been recognized in the past that 
both approximations should influence results in opposite 
directions, leading to a partial cancellation of errors, it 
is clear that the degree of cancellation depends on the 
thermodynamic conditions. Hence reliance on this can- 
cellation leads to predictions of limited validity. 
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